Environmental-variability-functional-traits
This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
# Load libraries
library(readr)
library(ggplot2)
library(nlme)
library(plyr)
library(gplots)##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowesslibrary(gstat)## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoolibrary(sp)
require(MuMIn)## Loading required package: MuMIn# SE function
se<-function(x) sqrt(var(x, na.rm=T)/length(x))
ggplotRegression <- function (fit) {
require(ggplot2)
ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) +
geom_point() +
stat_smooth(method = "lm", col = "red") +
labs(title = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5),
"Intercept =",signif(fit$coef[[1]],5 ),
" Slope =",signif(fit$coef[[2]], 5),
" P =",signif(summary(fit)$coef[2,4], 5)))
}Sceloporus_specimens = read_delim("Sceloporus_specimens.txt",
"\t", escape_double = FALSE, trim_ws = TRUE)## Parsed with column specification:
## cols(
## .default = col_double(),
## SeriesID = col_character(),
## Genus = col_character(),
## Species = col_character(),
## Date = col_character(),
## Area = col_character(),
## Site = col_character(),
## Time = col_time(format = ""),
## Sex = col_character(),
## Age10g = col_character(),
## VerbatimLocation = col_character()
## )## See spec(...) for full column specifications.# Select lizards >= 10 g only
Sceloporus_specimens_Age10g = Sceloporus_specimens[which(Sceloporus_specimens$Age10g=='Adult'),]
getwd()## [1] "/Users/vanw/Documents/Projects/Sceloporus projects/Morph evn MS/Data archive"# Descriptive stats
# SVL
summary(Sceloporus_specimens_Age10g$SVL_mm) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 64.00 71.75 76.00 76.66 81.00 93.00# Mass_g
summary(Sceloporus_specimens_Age10g$Mass_g)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.10 13.50 16.20 17.06 19.02 32.10# Dorsal scale counts
summary(Sceloporus_specimens_Age10g$Scales_dorsal)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 36.00 39.50 41.00 41.46 44.00 48.00 57# Dorsal color
summary(Sceloporus_specimens_Age10g$Color_Dorsal)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 18.98 34.78 37.75 37.52 41.21 48.78 7# data summary function
data_summary <- function(data, varname, groupnames){
require(plyr)
summary_func <- function(x, col){
c(mean = mean(x[[col]], na.rm=TRUE),
sd = sd(x[[col]], na.rm=TRUE))
}
data_sum<-ddply(data, groupnames, .fun=summary_func,
varname)
data_sum <- rename(data_sum, c("mean" = varname))
return(data_sum)
}# Table 1
df3.SVL <- data_summary(Sceloporus_specimens_Age10g, varname="SVL_mm",
groupnames=c("Site"))
write.table(df3.SVL, "df3.SVL.txt", sep="\t")
df3.dorsal.scales <- data_summary(Sceloporus_specimens_Age10g, varname="Scales_dorsal",
groupnames=c("Site"))
write.table(df3.dorsal.scales, "df3.dorsal.scales.txt", sep="\t")
df3.dorsal.color <- data_summary(Sceloporus_specimens_Age10g, varname="Color_Dorsal",
groupnames=c("Site"))
write.table(df3.dorsal.color, "df3.dorsal.color.txt", sep="\t")# SVL
par(mfrow=c(1,1))
hist(Sceloporus_specimens_Age10g$SVL_mm, breaks = 20)
# Summer
# Mean Temp of Warmest Quarter
par(mfrow=c(1,1))
par(mar=c(5,5,2,2))
library(lme4)## Loading required package: Matrix## Registered S3 method overwritten by 'lme4':
## method from
## predict.merMod MuMIn##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepreg1d.lmer = lmer(SVL_mm ~ BIO10_meanTwarmestQ + (1|BIO10_meanTwarmestQ), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1d.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO10_meanTwarmestQ + (1 | BIO10_meanTwarmestQ)
## Data: Sceloporus_specimens_Age10g
##
## AIC BIC logLik deviance df.resid
## 1073.6 1086.0 -532.8 1065.6 160
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0984 -0.5958 -0.1051 0.5966 2.5261
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIO10_meanTwarmestQ (Intercept) 17.47 4.180
## Residual 30.91 5.559
## Number of obs: 164, groups: BIO10_meanTwarmestQ, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 74.4157 5.1266 31.5532 14.52 1.63e-15 ***
## BIO10_meanTwarmestQ 0.1550 0.2673 31.6970 0.58 0.566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BIO10_mnTwQ -0.984anova(reg1d.lmer)plot(reg1d.lmer)
r.squaredGLMM(reg1d.lmer)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.## R2m R2c
## [1,] 0.005909285 0.3648981library(ggplot2)
reg1d.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO10_meanTwarmestQ, y = SVL_mm)) +
geom_point()
reg1d.1 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter (°C)", y = "Snout-Vent Length (mm)")
# Max Temp Warmest Month
library(lme4)
library(lmerTest)
reg1e.lmer = lmer(SVL_mm ~ BIO5_maxTwarmestMo + (1|BIO5_maxTwarmestMo), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1e.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO5_maxTwarmestMo + (1 | BIO5_maxTwarmestMo)
## Data: Sceloporus_specimens_Age10g
##
## AIC BIC logLik deviance df.resid
## 1071.4 1083.8 -531.7 1063.4 160
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2225 -0.6141 -0.1145 0.6227 2.4730
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIO5_maxTwarmestMo (Intercept) 20.66 4.545
## Residual 29.39 5.421
## Number of obs: 164, groups: BIO5_maxTwarmestMo, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 74.1033 6.6428 30.1156 11.155 3.21e-12 ***
## BIO5_maxTwarmestMo 0.1160 0.2228 30.5010 0.521 0.606
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BIO5_mxTwrM -0.990anova(reg1e.lmer)plot(reg1e.lmer)
r.squaredGLMM(reg1e.lmer)## R2m R2c
## [1,] 0.004640269 0.4155082library(ggplot2)
reg1d.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO5_maxTwarmestMo, y = SVL_mm)) +
geom_point()
reg1d.1 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month (°C)", y = "Snout-Vent Length (mm)")
# Mean Temp of Coldest Quarter
library(lme4)
library(lmerTest)
reg1b.lmer = lmer(SVL_mm ~ BIO11_meanTcoldestQ + (1|BIO11_meanTcoldestQ), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1b.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO11_meanTcoldestQ + (1 | BIO11_meanTcoldestQ)
## Data: Sceloporus_specimens_Age10g
##
## AIC BIC logLik deviance df.resid
## 1073.9 1086.3 -532.9 1065.9 160
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.46586 -0.63245 -0.03986 0.62416 3.02117
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIO11_meanTcoldestQ (Intercept) 17.38 4.168
## Residual 31.27 5.592
## Number of obs: 164, groups: BIO11_meanTcoldestQ, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 78.3629 1.2454 29.2064 62.921 <2e-16 ***
## BIO11_meanTcoldestQ -0.4895 0.3630 28.9627 -1.348 0.188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BIO11_mnTcQ -0.655anova(reg1b.lmer)plot(reg1b.lmer)
r.squaredGLMM(reg1b.lmer)## R2m R2c
## [1,] 0.03352773 0.3787038library(ggplot2)
reg1d.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO11_meanTcoldestQ, y = SVL_mm)) +
geom_point()
reg1d.1 + theme_classic() + labs(x = "Mean Temperature of Coldest Quarter (°C)", y = "Snout-Vent Length (mm)")
# Min temp coldest month
library(lme4)
library(lmerTest)
reg1c.lmer = lmer(SVL_mm ~ BIO6_minTcoldestMo + (1|BIO6_minTcoldestMo), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1c.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO6_minTcoldestMo + (1 | BIO6_minTcoldestMo)
## Data: Sceloporus_specimens_Age10g
##
## AIC BIC logLik deviance df.resid
## 1066.3 1078.7 -529.2 1058.3 160
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.92480 -0.63242 -0.08504 0.60222 2.45331
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIO6_minTcoldestMo (Intercept) 13.48 3.671
## Residual 30.49 5.521
## Number of obs: 164, groups: BIO6_minTcoldestMo, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 73.1123 1.8338 29.2611 39.87 <2e-16 ***
## BIO6_minTcoldestMo -0.8527 0.3452 28.7928 -2.47 0.0197 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BIO6_mnTclM 0.887anova(reg1c.lmer)plot(reg1c.lmer)
r.squaredGLMM(reg1c.lmer)## R2m R2c
## [1,] 0.09776281 0.3743755library(ggplot2)
reg1c.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO6_minTcoldestMo, y = SVL_mm)) +
geom_point() +
stat_smooth(method = "lm", col = "black")
reg1c.1 + theme_classic() + labs(x = "Minimum Temperature of Coldest Month (°C)", y = "Snout-Vent Length (mm)")
reg1c <- lm(Sceloporus_specimens_Age10g$SVL_mm~Sceloporus_specimens_Age10g$BIO6_minTcoldestMo)
summary(reg1c)##
## Call:
## lm(formula = Sceloporus_specimens_Age10g$SVL_mm ~ Sceloporus_specimens_Age10g$BIO6_minTcoldestMo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.8321 -5.3122 -0.2677 4.7100 15.8614
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 72.6509 1.1003 66.030
## Sceloporus_specimens_Age10g$BIO6_minTcoldestMo -0.8467 0.2044 -4.143
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Sceloporus_specimens_Age10g$BIO6_minTcoldestMo 5.51e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.679 on 162 degrees of freedom
## Multiple R-squared: 0.0958, Adjusted R-squared: 0.09022
## F-statistic: 17.16 on 1 and 162 DF, p-value: 5.507e-05Spatially-explicit models
note: some values my differ slightly from publication values due to variation in random jitter function of gps points. seed set for all spatially explict models in this example for reproducibility.
# Load dataset
library(readr)
library(ggplot2)
library(nlme)
library(plyr)
library(gplots)
# SE function
se<-function(x) sqrt(var(x, na.rm=T)/length(x))Sceloporus_VertNet = read_delim("Sceloporus_occidentalis_VertNet.txt",
"\t", escape_double = FALSE, trim_ws = TRUE)## Parsed with column specification:
## cols(
## SVL_mm = col_double(),
## decimallatitude = col_double(),
## decimallongitude = col_double(),
## SeriesID = col_character(),
## bio_5 = col_double(),
## bio_6 = col_double(),
## bio_10 = col_double(),
## bio_11 = col_double()
## )summary(Sceloporus_VertNet)## SVL_mm decimallatitude decimallongitude SeriesID
## Min. : 24.0 Min. :32.89 Min. :-122.7 Length:752
## 1st Qu.: 64.0 1st Qu.:37.35 1st Qu.:-119.7 Class :character
## Median : 73.0 Median :37.75 Median :-119.5 Mode :character
## Mean : 70.6 Mean :37.65 Mean :-119.3
## 3rd Qu.: 79.0 3rd Qu.:37.93 3rd Qu.:-118.6
## Max. :109.0 Max. :47.59 Max. :-115.4
## bio_5 bio_6 bio_10 bio_11
## Min. :198.0 Min. :-104.00 Min. :110.0 Min. :-29.00
## 1st Qu.:241.0 1st Qu.: -73.00 1st Qu.:144.0 1st Qu.: -9.00
## Median :285.5 Median : -54.00 Median :181.0 Median : 15.00
## Mean :277.5 Mean : -47.67 Mean :175.7 Mean : 17.71
## 3rd Qu.:307.0 3rd Qu.: -23.00 3rd Qu.:200.0 3rd Qu.: 37.00
## Max. :364.0 Max. : 58.00 Max. :260.0 Max. :122.00# Adjust temperature scale
Sceloporus_VertNet$bio_5 = Sceloporus_VertNet$bio_5/10
Sceloporus_VertNet$bio_6 = Sceloporus_VertNet$bio_6/10
Sceloporus_VertNet$bio_10 = Sceloporus_VertNet$bio_10/10
Sceloporus_VertNet$bio_11 = Sceloporus_VertNet$bio_11/10# Select lizards 64 g and larger only
Sceloporus_VertNet_65mm = Sceloporus_VertNet[which(Sceloporus_VertNet$SVL_mm>63.5),]summary(Sceloporus_VertNet_65mm)## SVL_mm decimallatitude decimallongitude SeriesID
## Min. : 64.00 Min. :33.68 Min. :-122.4 Length:577
## 1st Qu.: 71.00 1st Qu.:37.32 1st Qu.:-119.7 Class :character
## Median : 76.00 Median :37.74 Median :-119.5 Mode :character
## Mean : 76.47 Mean :37.64 Mean :-119.2
## 3rd Qu.: 81.00 3rd Qu.:37.93 3rd Qu.:-118.6
## Max. :109.00 Max. :47.59 Max. :-115.4
## bio_5 bio_6 bio_10 bio_11
## Min. :19.80 Min. :-10.400 Min. :11.00 Min. :-2.900
## 1st Qu.:23.70 1st Qu.: -7.300 1st Qu.:14.30 1st Qu.:-1.000
## Median :27.80 Median : -5.800 Median :17.50 Median : 1.200
## Mean :27.66 Mean : -5.035 Mean :17.46 Mean : 1.555
## 3rd Qu.:31.00 3rd Qu.: -2.800 3rd Qu.:20.00 3rd Qu.: 3.600
## Max. :36.40 Max. : 5.800 Max. :26.00 Max. :12.200hist(Sceloporus_VertNet_65mm$SVL_mm)
Mean Temperature of Warmest Quarter (bio10)
library(rgdal)## rgdal: version: 1.4-6, (SVN revision 841)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))
library(geoR)## Warning in system2("/usr/bin/otool", c("-L", shQuote(DSO)), stdout = TRUE):
## running command ''/usr/bin/otool' -L '/Library/Frameworks/R.framework/
## Resources/library/tcltk/libs//tcltk.so'' had status 1## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originpar(mfrow=c(1,1))
plot(v1)
library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
resamp=10, quiet = FALSE)## 1 of 10
2 of 10
3 of 10
4 of 10
5 of 10
6 of 10
7 of 10
8 of 10
9 of 10
10 of 10 plot(ncf.scor)
mod.nocorr = lm(SVL_mm ~ bio_10, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originplot(v1)
# jitter GPS points
set.seed(9999)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_10,
data = Sceloporus_VertNet_65mm,
correlation = corExp(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corExp_1a = SVL_mm.corExp
summary(SVL_mm.corExp_1a) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3733.249 3755.021 -1861.624
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 52.0663289 0.5255908
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 80.68188 3.406257 23.686374 0.0000
## bio_10 -0.24457 0.154141 -1.586665 0.1131
##
## Correlation:
## (Intr)
## bio_10 -0.92
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66309328 -0.68428635 -0.07113418 0.56200949 4.28756934
##
## Residual standard error: 7.751242
## Degrees of freedom: 577 total; 575 residual#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_1a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_10,
data = Sceloporus_VertNet_65mm,
correlation = corGaus(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corGaus_1b = SVL_mm.corGaus
summary(SVL_mm.corGaus_1b) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3738.805 3760.577 -1864.402
##
## Correlation Structure: Gaussian spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 31.4080885 0.5711817
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 81.69492 2.9314082 27.868831 0.0000
## bio_10 -0.29681 0.1365493 -2.173631 0.0301
##
## Correlation:
## (Intr)
## bio_10 -0.937
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.72288021 -0.69083161 -0.09976393 0.57325596 4.37452618
##
## Residual standard error: 7.605604
## Degrees of freedom: 577 total; 575 residualSVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_1b, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corGaus.var)
# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_10,
data = Sceloporus_VertNet_65mm,
correlation = corLin(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corLin_1c = SVL_mm.corLin
summary(SVL_mm.corLin_1c) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3734.05 3755.822 -1862.025
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.5094906 0.5646552
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 79.58240 3.348365 23.767543 0.0000
## bio_10 -0.18927 0.161908 -1.169009 0.2429
##
## Correlation:
## (Intr)
## bio_10 -0.963
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.70130731 -0.70563209 -0.04864043 0.61073152 4.46851351
##
## Residual standard error: 7.434682
## Degrees of freedom: 577 total; 575 residualSVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_1c, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corLin.var)
# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_10,
data = Sceloporus_VertNet_65mm,
correlation = corRatio(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corRatio_1d = SVL_mm.corRatio
summary(SVL_mm.corRatio_1d) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3738.134 3759.906 -1864.067
##
## Correlation Structure: Rational quadratic spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 32.7784921 0.5579583
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 81.58379 3.091487 26.389823 0.0000
## bio_10 -0.29214 0.138558 -2.108464 0.0354
##
## Correlation:
## (Intr)
## bio_10 -0.909
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.69715109 -0.68172480 -0.09456885 0.56885085 4.32730113
##
## Residual standard error: 7.692627
## Degrees of freedom: 577 total; 575 residualSVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_1d, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corRatio.var)
# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_10,
data = Sceloporus_VertNet_65mm,
correlation = corSpher(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corSpher_1e = SVL_mm.corSpher
summary(SVL_mm.corSpher_1e) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3735.045 3756.817 -1862.523
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 90.8478140 0.5282386
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 81.36789 3.276044 24.837239 0.0000
## bio_10 -0.27635 0.152577 -1.811218 0.0706
##
## Correlation:
## (Intr)
## bio_10 -0.938
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6882349 -0.6930102 -0.1013650 0.5508959 4.2819509
##
## Residual standard error: 7.75038
## Degrees of freedom: 577 total; 575 residualSVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_1e, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corSpher.var)
# calculate AIC to select best model
AIC(SVL_mm.corExp_1a, SVL_mm.corGaus_1b, SVL_mm.corLin_1c, SVL_mm.corRatio_1d, SVL_mm.corSpher_1e)# SVL_mm.corExp is best
# SVL_mm.corSpher & SVL_mm.corLin are close# best model
summary(SVL_mm.corExp_1a)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3733.249 3755.021 -1861.624
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 52.0663289 0.5255908
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 80.68188 3.406257 23.686374 0.0000
## bio_10 -0.24457 0.154141 -1.586665 0.1131
##
## Correlation:
## (Intr)
## bio_10 -0.92
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66309328 -0.68428635 -0.07113418 0.56200949 4.28756934
##
## Residual standard error: 7.751242
## Degrees of freedom: 577 total; 575 residual# other good models
summary(SVL_mm.corLin_1c)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3734.05 3755.822 -1862.025
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.5094906 0.5646552
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 79.58240 3.348365 23.767543 0.0000
## bio_10 -0.18927 0.161908 -1.169009 0.2429
##
## Correlation:
## (Intr)
## bio_10 -0.963
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.70130731 -0.70563209 -0.04864043 0.61073152 4.46851351
##
## Residual standard error: 7.434682
## Degrees of freedom: 577 total; 575 residualsummary(SVL_mm.corSpher_1e)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_10
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3735.045 3756.817 -1862.523
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 90.8478140 0.5282386
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 81.36789 3.276044 24.837239 0.0000
## bio_10 -0.27635 0.152577 -1.811218 0.0706
##
## Correlation:
## (Intr)
## bio_10 -0.938
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6882349 -0.6930102 -0.1013650 0.5508959 4.2819509
##
## Residual standard error: 7.75038
## Degrees of freedom: 577 total; 575 residual# all give same interpretationlibrary(ggplot2)
plot3 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_10, y = SVL_mm)) +
geom_point()
plot3 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter", y = "Snout-Vent Length (mm)")
reg3 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_10)
summary(reg3)##
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_10)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.983 -5.341 -0.544 4.345 33.770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.6603 1.4645 57.809 < 2e-16 ***
## Sceloporus_VertNet_65mm$bio_10 -0.4692 0.0821 -5.714 1.77e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.211 on 575 degrees of freedom
## Multiple R-squared: 0.05374, Adjusted R-squared: 0.05209
## F-statistic: 32.66 on 1 and 575 DF, p-value: 1.768e-08#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_1a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
#plot
library(ggplot2)
plot3.1 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_10, y = SVL_mm)) +
geom_point()
plot3.1 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter", y = "Snout-Vent Length (mm)")
Max Temp of Warmest Month (bio_5)
library(rgdal)
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))
library(geoR)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originpar(mfrow=c(1,1))
plot(v1)
library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
resamp=10, quiet = FALSE)## 1 of 10
2 of 10
3 of 10
4 of 10
5 of 10
6 of 10
7 of 10
8 of 10
9 of 10
10 of 10 plot(ncf.scor)
mod.nocorr = lm(SVL_mm ~ bio_5, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originplot(v1)
# jitter GPS points
set.seed(8888)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_5,
data = Sceloporus_VertNet_65mm,
correlation = corExp(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corExp_2a = SVL_mm.corExp
summary(SVL_mm.corExp_2a) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3734.34 3756.112 -1862.17
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 55.0865638 0.5194307
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 81.57497 4.241857 19.230957 0.0000
## bio_5 -0.19039 0.129958 -1.465044 0.1435
##
## Correlation:
## (Intr)
## bio_5 -0.946
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6321856 -0.6834408 -0.0421755 0.5776624 4.2676357
##
## Residual standard error: 7.804842
## Degrees of freedom: 577 total; 575 residual#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_2a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_5,
data = Sceloporus_VertNet_65mm,
correlation = corGaus(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corGaus_2b = SVL_mm.corGaus
summary(SVL_mm.corGaus_2b) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3740.048 3761.82 -1865.024
##
## Correlation Structure: Gaussian spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 31.9129845 0.5650735
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 82.77308 3.697749 22.384722 0.0000
## bio_5 -0.22993 0.115757 -1.986289 0.0475
##
## Correlation:
## (Intr)
## bio_5 -0.96
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.69063367 -0.69890215 -0.07982741 0.57979998 4.35722315
##
## Residual standard error: 7.649751
## Degrees of freedom: 577 total; 575 residualSVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_2b, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corGaus.var)
# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_5,
data = Sceloporus_VertNet_65mm,
correlation = corLin(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corLin_2c = SVL_mm.corLin
summary(SVL_mm.corLin_2c) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3735.683 3757.455 -1862.841
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.7132778 0.5630984
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 79.40055 4.317961 18.388433 0.0000
## bio_5 -0.11801 0.138808 -0.850169 0.3956
##
## Correlation:
## (Intr)
## bio_5 -0.978
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66430333 -0.70442457 -0.04417102 0.63742255 4.46110240
##
## Residual standard error: 7.452414
## Degrees of freedom: 577 total; 575 residualSVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_2c, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corLin.var)
# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_5,
data = Sceloporus_VertNet_65mm,
correlation = corRatio(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corRatio_2d = SVL_mm.corRatio
summary(SVL_mm.corRatio_2d) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3739.165 3760.937 -1864.583
##
## Correlation Structure: Rational quadratic spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.5808965 0.5514612
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 82.72663 3.833582 21.579463 0.00
## bio_5 -0.22947 0.116857 -1.963699 0.05
##
## Correlation:
## (Intr)
## bio_5 -0.939
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66625538 -0.68658515 -0.07497444 0.57688858 4.31019769
##
## Residual standard error: 7.740728
## Degrees of freedom: 577 total; 575 residualSVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_2d, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corRatio.var)
# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_5,
data = Sceloporus_VertNet_65mm,
correlation = corSpher(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corSpher_2e = SVL_mm.corSpher
summary(SVL_mm.corSpher_2e) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3736.25 3758.022 -1863.125
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 91.7092542 0.5251655
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 82.31749 4.146448 19.852530 0.0000
## bio_5 -0.21249 0.129907 -1.635743 0.1024
##
## Correlation:
## (Intr)
## bio_5 -0.961
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66136760 -0.70130724 -0.07443573 0.56013969 4.27534825
##
## Residual standard error: 7.776811
## Degrees of freedom: 577 total; 575 residualSVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_2e, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corSpher.var)
# calculate AIC to select best model
AIC(SVL_mm.corExp_2a, SVL_mm.corGaus_2b, SVL_mm.corLin_2c, SVL_mm.corRatio_2d, SVL_mm.corSpher_2e)# SVL_mm.corExp# best model
summary(SVL_mm.corExp_2a)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3734.34 3756.112 -1862.17
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 55.0865638 0.5194307
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 81.57497 4.241857 19.230957 0.0000
## bio_5 -0.19039 0.129958 -1.465044 0.1435
##
## Correlation:
## (Intr)
## bio_5 -0.946
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6321856 -0.6834408 -0.0421755 0.5776624 4.2676357
##
## Residual standard error: 7.804842
## Degrees of freedom: 577 total; 575 residual# other models with support
summary(SVL_mm.corLin_2c)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3735.683 3757.455 -1862.841
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.7132778 0.5630984
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 79.40055 4.317961 18.388433 0.0000
## bio_5 -0.11801 0.138808 -0.850169 0.3956
##
## Correlation:
## (Intr)
## bio_5 -0.978
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66430333 -0.70442457 -0.04417102 0.63742255 4.46110240
##
## Residual standard error: 7.452414
## Degrees of freedom: 577 total; 575 residualsummary(SVL_mm.corSpher_2e)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_5
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3736.25 3758.022 -1863.125
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 91.7092542 0.5251655
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 82.31749 4.146448 19.852530 0.0000
## bio_5 -0.21249 0.129907 -1.635743 0.1024
##
## Correlation:
## (Intr)
## bio_5 -0.961
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.66136760 -0.70130724 -0.07443573 0.56013969 4.27534825
##
## Residual standard error: 7.776811
## Degrees of freedom: 577 total; 575 residual# same interpretation#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_2a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
library(ggplot2)
plot4 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_5, y = SVL_mm)) +
geom_point()
plot4 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month", y = "Snout-Vent Length (mm)")
reg4 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_5)
summary(reg4)##
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.715 -5.440 -0.655 4.576 33.709
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.54793 1.97805 43.754 < 2e-16 ***
## Sceloporus_VertNet_65mm$bio_5 -0.36432 0.07066 -5.156 3.48e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.247 on 575 degrees of freedom
## Multiple R-squared: 0.04418, Adjusted R-squared: 0.04252
## F-statistic: 26.58 on 1 and 575 DF, p-value: 3.484e-07Winter
Mean Temp of Coldest Quarter (bio_11)
library(rgdal)
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))
library(geoR)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originpar(mfrow=c(1,1))
plot(v1)
library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
resamp=10, quiet = FALSE)## 1 of 10
2 of 10
3 of 10
4 of 10
5 of 10
6 of 10
7 of 10
8 of 10
9 of 10
10 of 10 plot(ncf.scor)
mod.nocorr = lm(SVL_mm ~ bio_11, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originplot(v1)
# jitter GPS points
set.seed(9999)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_11,
data = Sceloporus_VertNet_65mm,
correlation = corExp(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corExp_3a = SVL_mm.corExp
summary(SVL_mm.corExp_3a) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3729.706 3751.477 -1859.853
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 39.5944966 0.5678399
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.36239 1.3056903 59.25019 0.0000
## bio_11 -0.44154 0.1803243 -2.44860 0.0146
##
## Correlation:
## (Intr)
## bio_11 -0.498
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.92580993 -0.72430551 -0.09372455 0.56376069 4.45828937
##
## Residual standard error: 7.442991
## Degrees of freedom: 577 total; 575 residual#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_3a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_11,
data = Sceloporus_VertNet_65mm,
correlation = corGaus(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corGaus_3b = SVL_mm.corGaus
summary(SVL_mm.corGaus_3b) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3734.095 3755.867 -1862.048
##
## Correlation Structure: Gaussian spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 29.3955433 0.6049416
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.50703 1.1169734 69.39022 0.0000
## bio_11 -0.48784 0.1601221 -3.04665 0.0024
##
## Correlation:
## (Intr)
## bio_11 -0.516
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.97514173 -0.75594120 -0.09841634 0.55577894 4.49754792
##
## Residual standard error: 7.381887
## Degrees of freedom: 577 total; 575 residualSVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_3b, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corGaus.var)
# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_11,
data = Sceloporus_VertNet_65mm,
correlation = corLin(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corLin_3c = SVL_mm.corLin
summary(SVL_mm.corLin_3c) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3728.694 3750.466 -1859.347
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.4939287 0.5926691
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.41861 1.0533765 73.49567 0.0000
## bio_11 -0.45026 0.1707887 -2.63636 0.0086
##
## Correlation:
## (Intr)
## bio_11 -0.572
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.98416021 -0.74744628 -0.09743965 0.56886085 4.56579639
##
## Residual standard error: 7.262108
## Degrees of freedom: 577 total; 575 residualSVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_3c, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corLin.var)
# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_11,
data = Sceloporus_VertNet_65mm,
correlation = corRatio(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corRatio_3d = SVL_mm.corRatio
summary(SVL_mm.corRatio_3d) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3734.23 3756.001 -1862.115
##
## Correlation Structure: Rational quadratic spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 28.774088 0.594429
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.40568 1.3078337 59.18618 0.0000
## bio_11 -0.47775 0.1678583 -2.84615 0.0046
##
## Correlation:
## (Intr)
## bio_11 -0.459
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.94268900 -0.73327320 -0.08573498 0.56136307 4.47032918
##
## Residual standard error: 7.441609
## Degrees of freedom: 577 total; 575 residualSVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_3d, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corRatio.var)
# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_11,
data = Sceloporus_VertNet_65mm,
correlation = corSpher(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corSpher_3e = SVL_mm.corSpher
summary(SVL_mm.corSpher_3e) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3730.683 3752.455 -1860.342
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 53.5573820 0.5938406
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.46890 1.0923423 70.91998 0.000
## bio_11 -0.46464 0.1715471 -2.70851 0.007
##
## Correlation:
## (Intr)
## bio_11 -0.562
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.99123398 -0.75453630 -0.09968822 0.56182758 4.55617772
##
## Residual standard error: 7.277445
## Degrees of freedom: 577 total; 575 residualSVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_3e, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corSpher.var)
# calculate AIC to select best model
AIC(SVL_mm.corExp_3a, SVL_mm.corGaus_3b, SVL_mm.corLin_3c, SVL_mm.corRatio_3d, SVL_mm.corSpher_3e)#SVL_mm.corLinsummary(SVL_mm.corLin_3c)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3728.694 3750.466 -1859.347
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.4939287 0.5926691
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.41861 1.0533765 73.49567 0.0000
## bio_11 -0.45026 0.1707887 -2.63636 0.0086
##
## Correlation:
## (Intr)
## bio_11 -0.572
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.98416021 -0.74744628 -0.09743965 0.56886085 4.56579639
##
## Residual standard error: 7.262108
## Degrees of freedom: 577 total; 575 residual#best support# also look at these with approx. equal support
summary(SVL_mm.corExp_3a)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3729.706 3751.477 -1859.853
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 39.5944966 0.5678399
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.36239 1.3056903 59.25019 0.0000
## bio_11 -0.44154 0.1803243 -2.44860 0.0146
##
## Correlation:
## (Intr)
## bio_11 -0.498
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.92580993 -0.72430551 -0.09372455 0.56376069 4.45828937
##
## Residual standard error: 7.442991
## Degrees of freedom: 577 total; 575 residualsummary(SVL_mm.corSpher_3e)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_11
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3730.683 3752.455 -1860.342
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 53.5573820 0.5938406
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 77.46890 1.0923423 70.91998 0.000
## bio_11 -0.46464 0.1715471 -2.70851 0.007
##
## Correlation:
## (Intr)
## bio_11 -0.562
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.99123398 -0.75453630 -0.09968822 0.56182758 4.55617772
##
## Residual standard error: 7.277445
## Degrees of freedom: 577 total; 575 residual# same resultreg1 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_11)
summary(reg1)##
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_11)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.273 -5.166 -0.433 4.237 33.983
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.63003 0.32863 236.225 < 2e-16 ***
## Sceloporus_VertNet_65mm$bio_11 -0.74662 0.09529 -7.835 2.28e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.046 on 575 degrees of freedom
## Multiple R-squared: 0.09646, Adjusted R-squared: 0.09489
## F-statistic: 61.39 on 1 and 575 DF, p-value: 2.276e-14# abline(reg1, col = "black")#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_3a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
#plot
library(ggplot2)
plot1 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_11, y = SVL_mm)) +
geom_point() +
stat_smooth(method = "lm", col = "black")
plot1 + theme_classic() + labs(x = "Mean Temperature of Coldest Quarter", y = "Snout-vent Length (mm)")
Min Temp of Coldest Month (bio_6)
library(rgdal)
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))
library(geoR)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originpar(mfrow=c(1,1))
plot(v1)
library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
resamp=10, quiet = FALSE)## 1 of 10
2 of 10
3 of 10
4 of 10
5 of 10
6 of 10
7 of 10
8 of 10
9 of 10
10 of 10 plot(ncf.scor)
mod.nocorr = lm(SVL_mm ~ bio_6, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the originplot(v1)
# jitter GPS points
set.seed(9999)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_6,
data = Sceloporus_VertNet_65mm,
correlation = corExp(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corExp_4a = SVL_mm.corExp
summary(SVL_mm.corExp_4a) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3728.855 3750.627 -1859.428
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 35.0221945 0.5832795
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.29093 1.2020995 61.80098 0.0000
## bio_6 -0.48914 0.1829269 -2.67394 0.0077
##
## Correlation:
## (Intr)
## bio_6 0.468
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.88889142 -0.75134844 -0.08037974 0.53208840 4.53643131
##
## Residual standard error: 7.338495
## Degrees of freedom: 577 total; 575 residual#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_4a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_6,
data = Sceloporus_VertNet_65mm,
correlation = corGaus(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corGaus_4b = SVL_mm.corGaus
summary(SVL_mm.corGaus_4b) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3733.052 3754.824 -1861.526
##
## Correlation Structure: Gaussian spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 28.6290649 0.6175309
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.12319 1.0599616 69.93007 0.0000
## bio_6 -0.53011 0.1636516 -3.23927 0.0013
##
## Correlation:
## (Intr)
## bio_6 0.476
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.91497565 -0.77036119 -0.08412568 0.54836605 4.56258535
##
## Residual standard error: 7.307147
## Degrees of freedom: 577 total; 575 residualSVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_4b, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corGaus.var)
# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_6,
data = Sceloporus_VertNet_65mm,
correlation = corLin(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corLin_4c = SVL_mm.corLin
summary(SVL_mm.corLin_4c) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3727.277 3749.049 -1858.638
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.4908766 0.6021569
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.25792 1.0075366 73.70246 0.0000
## bio_6 -0.50057 0.1709142 -2.92878 0.0035
##
## Correlation:
## (Intr)
## bio_6 0.534
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.92996024 -0.76428663 -0.08470616 0.54689828 4.61822848
##
## Residual standard error: 7.208484
## Degrees of freedom: 577 total; 575 residualSVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_4c, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corLin.var)
# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_6,
data = Sceloporus_VertNet_65mm,
correlation = corRatio(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corRatio_4d = SVL_mm.corRatio
summary(SVL_mm.corRatio_4d) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3733.367 3755.139 -1861.683
##
## Correlation Structure: Rational quadratic spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 26.6737461 0.6093649
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.09573 1.2292637 60.27652 0.0000
## bio_6 -0.52290 0.1729929 -3.02267 0.0026
##
## Correlation:
## (Intr)
## bio_6 0.436
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.89383782 -0.75327970 -0.07532492 0.55633844 4.54478478
##
## Residual standard error: 7.34641
## Degrees of freedom: 577 total; 575 residualSVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_4d, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corRatio.var)
# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_6,
data = Sceloporus_VertNet_65mm,
correlation = corSpher(form = ~j.easting + j.northing,
nugget = TRUE,
value = c(50, 0.5)))
SVL_mm.corSpher_4e = SVL_mm.corSpher
summary(SVL_mm.corSpher_4e) ## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3729.347 3751.119 -1859.674
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 47.940820 0.609589
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.22969 1.0122070 73.33450 0.0000
## bio_6 -0.51405 0.1707524 -3.01051 0.0027
##
## Correlation:
## (Intr)
## bio_6 0.529
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.9480198 -0.7768292 -0.0899644 0.5539295 4.6365309
##
## Residual standard error: 7.177685
## Degrees of freedom: 577 total; 575 residualSVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_4e, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corSpher.var)
# calculate AIC to select best model
AIC(SVL_mm.corExp_4a, SVL_mm.corGaus_4b, SVL_mm.corLin_4c, SVL_mm.corRatio_4d, SVL_mm.corSpher_4e)# SVL_mm.corLin# best support
summary(SVL_mm.corLin_4c)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3727.277 3749.049 -1858.638
##
## Correlation Structure: Linear spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 33.4908766 0.6021569
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.25792 1.0075366 73.70246 0.0000
## bio_6 -0.50057 0.1709142 -2.92878 0.0035
##
## Correlation:
## (Intr)
## bio_6 0.534
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.92996024 -0.76428663 -0.08470616 0.54689828 4.61822848
##
## Residual standard error: 7.208484
## Degrees of freedom: 577 total; 575 residual# also good
summary(SVL_mm.corExp_4a)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3728.855 3750.627 -1859.428
##
## Correlation Structure: Exponential spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 35.0221945 0.5832795
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.29093 1.2020995 61.80098 0.0000
## bio_6 -0.48914 0.1829269 -2.67394 0.0077
##
## Correlation:
## (Intr)
## bio_6 0.468
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.88889142 -0.75134844 -0.08037974 0.53208840 4.53643131
##
## Residual standard error: 7.338495
## Degrees of freedom: 577 total; 575 residualsummary(SVL_mm.corSpher_4e)## Generalized least squares fit by REML
## Model: SVL_mm ~ bio_6
## Data: Sceloporus_VertNet_65mm
## AIC BIC logLik
## 3729.347 3751.119 -1859.674
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~j.easting + j.northing
## Parameter estimate(s):
## range nugget
## 47.940820 0.609589
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 74.22969 1.0122070 73.33450 0.0000
## bio_6 -0.51405 0.1707524 -3.01051 0.0027
##
## Correlation:
## (Intr)
## bio_6 0.529
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.9480198 -0.7768292 -0.0899644 0.5539295 4.6365309
##
## Residual standard error: 7.177685
## Degrees of freedom: 577 total; 575 residual# same resultreg2 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_6)
summary(reg2)##
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.457 -5.200 -0.272 4.309 34.320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 72.24936 0.56118 128.746 <2e-16 ***
## Sceloporus_VertNet_65mm$bio_6 -0.83810 0.09544 -8.781 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.961 on 575 degrees of freedom
## Multiple R-squared: 0.1182, Adjusted R-squared: 0.1167
## F-statistic: 77.11 on 1 and 575 DF, p-value: < 2.2e-16#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_4a, type = "normalized"))## variog: computing omnidirectional variogramplot(SVL_mm.corExp.var)
# plot
library(ggplot2)
plot2 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_6, y = SVL_mm)) +
geom_point() +
stat_smooth(method = "lm", col = "black")
plot2 + theme_classic() + labs(x = "Minimum Temperature of Coldest Month", y = "Snout-vent Length (mm)")
# Scales
# Load libraries
library(readr)
library(ggplot2)
library(nlme)
library(plyr)
library(gplots)
library(gstat)
library(sp)
# SE function
se<-function(x) sqrt(var(x, na.rm=T)/length(x))Sceloporus_specimens = read_delim("Sceloporus_specimens.txt",
"\t", escape_double = FALSE, trim_ws = TRUE)## Parsed with column specification:
## cols(
## .default = col_double(),
## SeriesID = col_character(),
## Genus = col_character(),
## Species = col_character(),
## Date = col_character(),
## Area = col_character(),
## Site = col_character(),
## Time = col_time(format = ""),
## Sex = col_character(),
## Age10g = col_character(),
## VerbatimLocation = col_character()
## )## See spec(...) for full column specifications.summary(Sceloporus_specimens)## SeriesID LACM Genus Species
## Length:217 Min. :188168 Length:217 Length:217
## Class :character 1st Qu.:188222 Class :character Class :character
## Mode :character Median :188278 Mode :character Mode :character
## Mean :188278
## 3rd Qu.:188333
## Max. :188387
##
## Date Area Site
## Length:217 Length:217 Length:217
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## Time Sex Latitude Longitude
## Length:217 Length:217 Min. :34.10 Min. :-122.3
## Class1:hms Class :character 1st Qu.:34.37 1st Qu.:-118.6
## Class2:difftime Mode :character Median :37.26 Median :-118.4
## Mode :numeric Mean :36.65 Mean :-118.7
## 3rd Qu.:37.36 3rd Qu.:-117.7
## Max. :41.28 Max. :-116.9
##
## PosErr Elevation SVL_mm Age10g
## Min. :2.000 Min. : 973 Min. :42.00 Length:217
## 1st Qu.:2.000 1st Qu.:1239 1st Qu.:66.00 Class :character
## Median :2.000 Median :1532 Median :73.00 Mode :character
## Mean :2.451 Mean :1655 Mean :72.55
## 3rd Qu.:3.000 3rd Qu.:2216 3rd Qu.:79.00
## Max. :5.000 Max. :2456 Max. :93.00
## NA's :2
## Mass_g VerbatimLocation Color_Dorsal Scales_dorsal
## Min. : 2.9 Length:217 Min. :18.98 Min. :36.00
## 1st Qu.:10.3 Class :character 1st Qu.:33.51 1st Qu.:39.00
## Median :14.7 Mode :character Median :36.75 Median :41.00
## Mean :14.8 Mean :36.80 Mean :41.42
## 3rd Qu.:17.7 3rd Qu.:40.49 3rd Qu.:43.00
## Max. :32.1 Max. :48.78 Max. :48.00
## NA's :8 NA's :76
## BIO12_annualPrecip BIO5_maxTwarmestMo BIO6_minTcoldestMo Aridity_index
## Min. : 161.0 Min. :22.4 Min. :-8.500 Min. : 13.08
## 1st Qu.: 365.0 1st Qu.:25.0 1st Qu.:-6.500 1st Qu.: 47.97
## Median : 548.0 Median :28.6 Median :-4.800 Median : 72.21
## Mean : 628.4 Mean :28.9 Mean :-4.573 Mean : 87.11
## 3rd Qu.: 743.0 3rd Qu.:31.6 3rd Qu.:-2.100 3rd Qu.:116.85
## Max. :1760.0 Max. :35.7 Max. : 0.400 Max. :242.62
##
## BIO10_meanTwarmestQ BIO11_meanTcoldestQ
## Min. :13.20 Min. :-2.200
## 1st Qu.:16.00 1st Qu.: 1.200
## Median :18.10 Median : 2.300
## Mean :18.57 Mean : 2.447
## 3rd Qu.:21.70 3rd Qu.: 3.900
## Max. :25.10 Max. : 7.500
## # Select lizards > 10 g only
Sceloporus_specimens_Age10g = Sceloporus_specimens[which(Sceloporus_specimens$Age10g=='Adult'),]# select lizards with scale data only
df.complete = subset(Sceloporus_specimens_Age10g, !(is.na(Scales_dorsal)))
plot(df.complete$Scales_dorsal, df.complete$SVL_mm)
SVL.Scales = lm(df.complete$SVL_mm~df.complete$Scales_dorsal)
abline(SVL.Scales, col = "red")
SVL.Scales.resid = resid(SVL.Scales)
SVL.Scales.predict = predict(SVL.Scales)
length(resid(SVL.Scales))## [1] 107length(predict(SVL.Scales))## [1] 107par(mfrow=c(1,1))
plot(predict(SVL.Scales), resid(SVL.Scales))
abline(0,0)
# POSITIVE residuals are larger scales, NEGATIVE resuduals are smaller scalesresid.Scales_dorsal = resid(SVL.Scales)
df.complete$resid.Scales_dorsal = resid.Scales_dorsal# Correlations among scale traits
df.complete$Scales_size = (df.complete$SVL_mm/df.complete$Scales_dorsal)
df.complete$Scale.count = df.complete$Scales_dorsal
df.complete$Scale.size = df.complete$Scales_size
# size vs number of scales
cor.test(df.complete$SVL_mm, df.complete$Scales_dorsal)##
## Pearson's product-moment correlation
##
## data: df.complete$SVL_mm and df.complete$Scales_dorsal
## t = 4.5891, df = 105, p-value = 1.238e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2372919 0.5554874
## sample estimates:
## cor
## 0.408736# size vs mean scale size
cor.test(df.complete$SVL_mm, df.complete$Scale.size)##
## Pearson's product-moment correlation
##
## data: df.complete$SVL_mm and df.complete$Scale.size
## t = 10.721, df = 105, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6178455 0.8026060
## sample estimates:
## cor
## 0.7229049SIMPLE CORRELATION FIGURES
library(ggplot2)
modelR1 = ggplot(df.complete, aes(x = SVL_mm, y = Scales_dorsal)) +
geom_point() +
stat_smooth(method = "lm", col = "black")
modelR1 + theme_classic() + labs(x = "Snout-Vent Length (mm)", y = "Number of Dorsal Scales")
library(ggplot2)
modelR2 = ggplot(df.complete, aes(x = SVL_mm, y = Scale.size)) +
geom_point() +
stat_smooth(method = "lm", col = "black")
modelR2 + theme_classic() + labs(x = "Snout-Vent Length (mm)", y = "Mean Dorsal Scale Size (mm)")
# Mean Temperature of Warmest Quarter
reg3d2.lmer.BIO10 = lmer(resid.Scales_dorsal ~ BIO10_meanTwarmestQ + (1|BIO10_meanTwarmestQ), REML = F, na.action=na.omit, data=df.complete)
#summary(reg3d2.lmer.BIO10)
anova(reg3d2.lmer.BIO10)plot(reg3d2.lmer.BIO10)
r.squaredGLMM(reg3d2.lmer.BIO10)## R2m R2c
## [1,] 0.003235912 0.1094919library(ggplot2)
plot2 = ggplot(df.complete, aes(x = BIO10_meanTwarmestQ, y = SVL_mm)) +
geom_point()
plot2 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter (°C)", y = "Snout-Vent Length (mm)")
# Maximum Temperature of Warmest Month
reg3d2.lmer.BIO5 = lmer(resid.Scales_dorsal ~ BIO5_maxTwarmestMo + (1|BIO5_maxTwarmestMo), REML = F, na.action=na.omit, data=df.complete)
#summary(reg3d2.lmer.BIO5)
anova(reg3d2.lmer.BIO5)plot(reg3d2.lmer.BIO5)
r.squaredGLMM(reg3d2.lmer.BIO5)## R2m R2c
## [1,] 0.003710512 0.1119964library(ggplot2)
plot2 = ggplot(df.complete, aes(x = BIO5_maxTwarmestMo, y = SVL_mm)) +
geom_point()
plot2 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month (°C)", y = "Snout-Vent Length (mm)")
# Aridity
reg3d2.lmer = lmer(resid.Scales_dorsal ~ Aridity_index + (1|Aridity_index), REML = F, na.action=na.omit, data=df.complete)
summary(reg3d2.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: resid.Scales_dorsal ~ Aridity_index + (1 | Aridity_index)
## Data: df.complete
##
## AIC BIC logLik deviance df.resid
## 686.9 697.5 -339.4 678.9 103
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06675 -0.58632 0.07179 0.56288 2.25282
##
## Random effects:
## Groups Name Variance Std.Dev.
## Aridity_index (Intercept) 1.489 1.220
## Residual 32.000 5.657
## Number of obs: 107, groups: Aridity_index, 30
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.044325 1.013973 22.133786 2.016 0.0561 .
## Aridity_index -0.024633 0.009427 23.243250 -2.613 0.0155 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Aridity_ndx -0.792anova(reg3d2.lmer)r.squaredGLMM(reg3d2.lmer)## R2m R2c
## [1,] 0.07119116 0.1124801library(ggplot2)
reg3d2 = ggplot(df.complete, aes(x = Aridity_index, y = resid.Scales_dorsal)) +
geom_point() +
stat_smooth(method = "lm", col = "black")
reg3d2 + theme_classic() + labs(x = "Aridity Index", y = "Dorsal Scale Size (resid.)")
# Dorsal color
par(mfrow=c(1,1))
hist(Sceloporus_specimens_Age10g$Color_Dorsal, breaks = 20)
# Summer
# Mean Temp of Warmest Quarter
par(mfrow=c(1,1))
par(mar=c(5,5,2,2))
library(lme4)
library(lmerTest)
reg1d2.lmer = lmer(Color_Dorsal ~ BIO10_meanTwarmestQ + (1|BIO10_meanTwarmestQ), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1d2.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: Color_Dorsal ~ BIO10_meanTwarmestQ + (1 | BIO10_meanTwarmestQ)
## Data: Sceloporus_specimens_Age10g
##
## AIC BIC logLik deviance df.resid
## 933.0 945.2 -462.5 925.0 153
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.63890 -0.47611 -0.01983 0.61106 2.09243
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIO10_meanTwarmestQ (Intercept) 10.86 3.296
## Residual 16.46 4.057
## Number of obs: 157, groups: BIO10_meanTwarmestQ, 33
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 29.0645 3.9639 28.2431 7.332 5.24e-08 ***
## BIO10_meanTwarmestQ 0.4148 0.2071 28.6045 2.003 0.0548 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BIO10_mnTwQ -0.984anova(reg1d2.lmer)plot(reg1d2.lmer)
r.squaredGLMM(reg1d2.lmer)## R2m R2c
## [1,] 0.06950942 0.4394815library(ggplot2)
reg1d2 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO10_meanTwarmestQ, y = Color_Dorsal)) +
geom_point()
reg1d2 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter (°C)", y = "Dorsal Lightness")## Warning: Removed 7 rows containing missing values (geom_point).
# Max Temp Warmest Month
library(ggplot2)
reg1e2 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO5_maxTwarmestMo, y = Color_Dorsal)) +
geom_point()
reg1e2 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month (°C)", y = "Dorsal Lightness")## Warning: Removed 7 rows containing missing values (geom_point).
library(lme4)
library(lmerTest)
reg1e2.lmer = lmer(Color_Dorsal ~ BIO5_maxTwarmestMo + (1|BIO5_maxTwarmestMo), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1e2.lmer)## Linear mixed model fit by maximum likelihood . t-tests use
## Satterthwaite's method [lmerModLmerTest]
## Formula: Color_Dorsal ~ BIO5_maxTwarmestMo + (1 | BIO5_maxTwarmestMo)
## Data: Sceloporus_specimens_Age10g
##
## AIC BIC logLik deviance df.resid
## 938.7 950.9 -465.3 930.7 153
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9847 -0.4566 -0.0250 0.5928 2.3138
##
## Random effects:
## Groups Name Variance Std.Dev.
## BIO5_maxTwarmestMo (Intercept) 10.42 3.228
## Residual 17.15 4.141
## Number of obs: 157, groups: BIO5_maxTwarmestMo, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 32.0813 4.8770 31.7622 6.578 2.13e-07 ***
## BIO5_maxTwarmestMo 0.1714 0.1641 32.5807 1.044 0.304
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## BIO5_mxTwrM -0.990anova(reg1e2.lmer)plot(reg1e2.lmer)
r.squaredGLMM(reg1e2.lmer)## R2m R2c
## [1,] 0.01758458 0.3888856